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Modeling the thermal evolution of 
enzyme-created bubbles in DNA 
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The formation of bubbles in nucleic acids (NAs) are fundamental in many biological processes 
such as DNA replication, recombination, telomeres formation, nucleotide excision repair, as 
well as RNA transcription and splicing. These precesses are carried out by assembled complexes 
with enzymes that separate selected regions of NAs. Within the frame of a nonlinear dynamics 
approach we model the structure of the DNA duplex by a nonlinear network of coupled 
oscillators. We show that in fact from certain local structural distortions there originate 
oscillating localized patterns, that is radial and torsional breathers, which are associated with 
localized H-bond deformations, being reminiscent of the replication bubble. We further study 
the temperature dependence of these oscillating bubbles. To this aim the underlying nonlinear 
oscillator network of the DNA duplex is brought in contact with a heat bath using the Nose- 
Hoover-method. Special attention is paid to the stability of the oscillating bubbles under 
the imposed thermal perturbations. It is demonstrated that the radial and torsional breathers, 
sustain the impact of thermal perturbations even at temperatures as high as room temperature. 
Generally, for nonzero temperature the H-bond breathers move coherently along the double 
chain whereas at T = standing radial and torsional breathers result. 

Keywords: DNA, enzymes, bubbles, breathers. PACS: 87.-15. v, 63.20.Kr, 63.20.Ry 



1. Introduction 

The vital life processes of nucleic acids (NAs) have 
become increasingly addressed and intensively studied 
over the last years by biologists, chemists, mathemati- 
cians and physicists. The interest of the latter is mainly 
focused on the actual dynamical processes involved in 
information retrieval from the NAs duplexes associated 
with structural transitions to single-stranded sequences. 
Consequently, the two complementary strands of double 
helical DNA must be separated (or at least locally 
melted) and rearranged in order to yield single-stranded 
sequences to be used by other molecules, requiring the 
breaking up of the hydrogen bonds in the corresponding 
region of the double helix. Therefore, the opening of 
a select region of the twisted (Watson- Crick) DNA 
double helix, begins with a partial unwinding at an 
area called the replication fork and is observed as a 
bubble. In general, the unwinding of NAs is achieved by 
a type of enzyme b elonging to the larg e family of nucleic 
acid hehcases fsee lStrver et al.ir2nn2|) (more than sixty 
different types). 

They are ubiquitous and versatile enzymes that, 
in conjunction with other components of the macro- 
molecular machines, carry out important biological 
processes such as DNA replication, DNA recombina- 
tion, nucleotide excision repair (i.e. UV damage), RNA 
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editing and splicing, transfer of single-stranded nucleic 
acid (ss NA) to other ssNA, protein o r release into 
solution ijvon Hippel fc Delagouttel200ll) and formation 
of T-loop telomeres (jWu fc Hicksonll200l|) . Nucleic acid 
helicases have several structural varieties (monomeric, 
dimeric, trimeric, tetrameric and closed hexameric), but 
all of them use the hydrolysis of nucleoside triphos- 
phate (NTP) to nucleoside diphosphate (NDP) as the 
preferred source of energy. Our present study deals 
with a nonlinear dynamics model of the formation of 
oscillating bubbles in regions of duplex DNA acted 
upon by an enzyme. Utilizing models based on nonlinear 
lattice dynamics to study the formation of DN A open- 
ing has been of great interest lately (Englander et alJ 
llQSOt lAgarwal fc Henni^l2003|) . The existence of local- 
ized modes, such as solitons and breathers (describ- 
ing energy localization and coherent transport), makes 
the nonlinear lattice approaches appealing. With 
view to vibrational excitation in DNA the P eyrard- 
Bishop CPB ) model llPevrard fc_Bishodll989l)_ and its 
successors ( Dauxois ct al. 1993; Barbi ct al. 1999a; 
Cocco fc Monasson, .1999: .Barbi et al., ■1999b: .Barbii 
1998|) have been successfully applied to describe moving 



localized excitations (breathers) which reproduce typi- 
cal features of the DNA opening dynamics such as the 
magnitude of the amplitudes and the time scale of the 
oscillating 'bubble' preceding full strand separation. 
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We are interested in the bubble formation process 
initiated by structural deformations of selected regions 
of the parental DNA duplex that serves as a template 
for the replication. The process begins when the repli- 
cation apparatus identifies the starting point and then 
gets combined to it. The replication apparatus is and 
assembled complex that forms the replication fork and 
opens it directionally. During this process the helical 
DNA is unwound and replicated. One of the components 
of that complex is the helicase, an enzyme unable to 
recognize the origin of the replication per se, and which 
requires the particip ation of specific proteins to lead it 
to the initiation site IJDelagoutte fc von HiT^De]ll200,1^ . 

It is assumed that this enzyme operates in the way, 
that a local unwinding in combination with (rather 
small) stretchings of the H-bonds in this region occurs. 
We aim to demonstrate that there are indeed initial 
deformations that give rise to localized vibrations, 
constrained to a region of DNA, matching the prop- 
erties of oscillating bubbles observed experimentally in 
DNA. These oscillating bubbles with their temporarily 
extended but yet unbroken H-bonds serve as the pre- 
cursors to the replication bubble. Furthermore, we are 
also interested whether the stable radial and torsional 
breathers persist under imposed thermal perturbations. 

2. Methods 

The nonlinear oscillator network model for the DNA 
double helix used in this paper is explained in 
detail in (Barbi et al. 1999a: Cocco & Monasson 1999; 



lAgarwal fc Hennie 2003: Hcnnig fc Archilla 2004) . The 
equilibrium position of each base within the duplex 
configuration is described in a Cartesian coordinate 
system by a;„|, y„ ^ and z^^-. The index pair (n, i) 
labels the n-th base on the i— th strand with i = 1, 2 
and 1 ^ ri ^ A'^ , where N is the number of base pairs 
considered. Displacements of the bases from their equi- 
librium positions are denoted by Xn,i, yn,i and z„_i. The 
potential energy taking into account the interactions 
between the bases consists of four parts. The potential 
energy of the hydrogen bond within a base pair is 
modeled typically by a Morse potential 



v;: = D„ 



exp 



/ a 
[-2' 



1 



(1) 



where the variables d„ describe dynamical deviations 
of the hydrogen b onds from their equilibri um lengths 
do (for details see iHennig fc Archillal l2004|) . The site- 
dependent depth of the Morse potential, £>„, depends 
on the number of involved hydrogen bonds for the two 
different pairings in DNA, namely the G-C and the A- 
T pairs. The former pair includes three hydrogen bonds 
while the latter includes only two. a~^ is a measure of 
the potential-well width. 

The energies of the rather strong and rigid covalent 
bonds between the nucleotides n and n — 1 on the i — th 
strand are modeled by harmonic potential terms 
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and ln,i describes the deviations from the equilibrium 
distance between two adjacent bases on the same 
strand. K is the elasticity coefficient. 



Effects of stacking, which impede that, due to 
the backbone rigidity, one base slides over another 
(ptrver ct al. 2002) are incorporated in the following 
potential terms 
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The supposedly small deformations in longitudinal 
direction can be modeled by harmonic elasticity poten- 
tial terms given by 

C 






The kinetic energy of a nucleotide is determined by 
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where ra is the mass and Vnf'^ 
(x, y, z)— component of the momentum. 
The model Hamiltonian reads then as 
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and the summation in (jSJ is performed over all 
nucleotides and the two strands. 

For a typical equilibrium configuration the rotation 
angle, by which each base is rotated around the central 
axis, is given by 6o = 36°, the distance between base 
pair planes is ft, = 3.4 A, and the inter-base distance 
(the diameter of the helix) is do = 20 A. For the aver- 
age mass of one nucleotide we use M = 300 amu = 
4.982 X lO^^^kg. The arbitrary base pair sequence is 
reflected in randomly distributed bi-valued H-bond 
couphng strengths Do (for A-T) and Di =2 Do (for G- 
C). Note that the ratio 1.5 is often used due to the 
fact that the G-C base pair has 3 hydrogen bridges 
and the A-T two of them, but the energy depends on 
the angles and distances of the H-bonds and on the 
polarity of the bases among other facto rs. Quantum 
chemical calculations l|SDoner et alJl200ll) lead to the 
ratio 2 as used here. However, the ratio 1.5 lead to 
qualitatively similar results. The aperiodic arrangement 
of the two different base-pairs, coding the genetic 
information, renders the DNA chain to be of A-B- 
disorder type. Following Barbi et al l)l999at) we set 
a = 4.45A-i, Do = OMeY, and K=lFeVA~^. The 
value of the parameter S i s reasonably taken as S" = 2K 
i|Cocco fc Monassor Il999^ and a credible val ue for C is 
given hv C = S/W i|Agarwal fc Hennig'200 3') . 

With the time scaled as t^ ^/Doo^ jm t one passes 
to a dimensionless formulation with quantities: 

S^) Jv) „(^) 
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Subsequently, the tildes are dropped. 

As for the temperature dependence we assume that 
the system is in contact with a heat bath which based 
on molecular dynamics t echniques is simulated wit h 
the Nose-Hoover method ()Noselll984albl: iHoovedllOSSJ) . 
According to this method the DNA lattice system is 
coupled to two additional variables, Ps and s, such that 
the equations read as 
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Note that the heat bath variable s does not enter 
explicitly the equations of motion regarding the actual 
dynamics of the bases. On the other hand H — H + 
p1/{2Q) + NksThi^s) is an augmented Hamiltonian 
the conservation of which can be used to assure accu- 
racy of the numerical computations. The parameter 
Q determines the time scale of the thermostat with 
temperature T and ks is the Boltzmann constant. 

The values of the scaled parameters are given hy K — 
0.683, ro= 44.50, /i = 15.13 and Zo = 31.39. One time 
unit of the scaled time corresponds to 0.198 ps of the 
physical time. 

We turn now to the study of the creation of bubbles 
in DNA. The starting point is a DNA molecule for 
which a certain segment experiences initially angular 
and radial deformations due to the action of some 
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Figure 1. The initial radial distortions of the double helix. 
The inter-base distance dn(0) in A. 
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Figure 2. The initial angular distortions of the double helix. 



enzyme to which a region of the DNA is bound. In order 
to simulate the deforming action of enzymes we assume 
that initially a number of consecutive sites in the center 
of the DNA lattice (hereafter referred to as the central 
region) are exerted to forces acting in angular and radial 
direction such that in this region the molecule experi- 
ences twist reduction together with radial stretchings. 
These structural deformations can be extended over 
a region encompassing up to thirty base pairs and as 
it is going to be demonstrated give rise to the for- 
mation of H-bridge breather solutions (extending over 
15 — 20 base pairs) reproducing the oscillating 'bubbles' 
observed for the DNA-opening process (Barbi 1998) . In 
Figs. Hand 13 the localized initial distortions are shown. 
Both the angular and radial deformation patterns are 
bell-shaped, due to the fact that the enzymatic force 
is exerted locally, i.e. in the extreme case to a single 
base pair only for which the H-bond is deformed ^ . The 
first one being of non-positive amplitudes is linked with 
reduced twist while the latter one with non-negative 
amplitudes is associated with radial stretchings. The 
radial and angular deformation patterns are centered at 
the central lattice site (base pair) at which the H-bond 
stretching and twist angle reduction is at maximum. At 
either side of the central site the amplitudes approach 
progressively zero. The deformation energy amounts to 
0.0362 eV. The set of coupled equations (IHIl-llIlIl is 
integrated with the help of a fourth-order Runge-Kutta 

-l-It turns out that this shape is also more stable in a ther mal bath 
than a rectangular shape as in lHennig fc Archillall20n4L 
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Figure 3. The case T — 0: Spatio-temporal pattern of the 
inter-base distance d„(f) in A. Parameters; C = 0.126, K = 
0.63, 5=1.26, Do = l, Di=2. 



method. For the simulation the DNA lattice consists of 
799 sites and open boundary conditions were imposed. 
The same initial conditions are used both at zero and 
nonzero temperatures. 

3. Results 

First of all we consider the zero temperature regime. 
In Fig. 13 we depict the spatio-temporal evolution 
of the distance, d„(t), between two bases of a base 
pair, which measures the variation of the length of 
the corresponding hydrogen bond expressed in A. The 
dynamics starts from a non-equilibrium configuration so 
that energy redistribution takes place. Most noticeably, 
this is established in the immediate emission of (small- 
amplitude) phonon waves from either side of the local- 
ized radial pattern situated around the central lattice 
site. These primary phonon waves travel uniformly 
towards the ends of the lattice with uniform velocity. In 
the approach of an equilibrium regime small amount of 
excitation energy is dispersed in the form of secondary 
phonon waves in the rest of the DNA lattice. However, 
the vast majority of the excitation energy remains 
contained in the central region. 

We observe that the localized bell-shaped radial 
pattern is preserved and performs periodic oscillations 
in time, viz. a radial breather is formed. This periodi- 
cally oscillating pattern of the radial variable, dn(t), is 
attributed to successive stretchings and compressings of 
the hydrogen bonds. 

Note that the stretching of a base pair distance is 
larger than the compressio n typical for the evo lution in 
a Morse potential (see also lBarbi et al.lll999al) . 

With regard to the associated pattern of the devia- 
tions of the twist angles from their equilibrium values, 
9 nit) — n0Q, expressed in rad, we find that in an initial 
phase the amplitudes increase. With the emanation of 
the primary radial phonon waves from the localized 
radial pattern (cf. Fig. ISJ there is linked a split up of 
the standing angular bell-shaped pattern so that two 
primary torsional waves of negative localized angular 
deformations, symmetrically to the central site, get 
produced which travel in unison with the radial phonon 
waves in the direction of the lattice ends. Since the 
radial phonon waves possess positive amplitudes corre- 
sponding to H-bond stretching the region of the double 
helix traversed by them experiences reductions of the 




Figure 4. The case T — Q: Spatio-temporal pattern of the 
twist angle expressed in rad. Parameters as in Fig. El and 
heat bath parameter Q — 10. 



twist angle, i.e. 9n(t) — uOq is negative. Similarly, the 
second radial phonon waves are connected with two 
small negative-amplitude angular waves following the 
primary ones away from the central region. Around the 
central lattice site the (initially negative) amplitudes of 
the twist angle grow steadily. Finally, after nearly 150 
time units the twist angles of the helix in the breather 
region become positive being equivalent to increased 
twist in a region comprising 20 lattice sites to either side 
of the central site. Notice that the process of amplitude 
breathing of the torsional part of the lattice proceeds 
with a longer period than the breathing of the radial 
localized pattern. This behavior of 9n{t) — n9o for zero 
temperature is shown in Fig. ^ 

Remarkably, in the non-zero temperature regime, 
when thermal energy is injected into the DNA lattice, 
the radial breather basically remains in localized shape. 
Moreover, its amplitudes and energy grow. Therefore, 
the breather extracts energy from the thermal bath. 
Simultaneously we observe that excess energy, not to 
be carried by the radial breather, is ejected from the 
central region where it disperses into wider parts of the 
lattice as phonons. Nevertheless, a radial breather of 
fairly large amplitude prevails over the small-amplitude 
noisy background caused by the heat bath on the DNA 
lattice. Interestingly, under the impact of thermal per- 
turbations the radial breather starts to move towards 
one end of the DNA lattice in coherent fashion with 
maximum amplitude being larger then the one in the 
T — case. This feature is illustrated for the case T = 
100 in Fig. Correspondingly, the propagating radial 
breather is accompanied by a torsional breather leav- 
ing the double helix consecutively in over-twisted and 
under-twisted shape. The direction of the movement of 
the breather is determined by the actual realization of 
the random sequence of the DNA base pairs. This is 
reflected in our model in bi-valued dissociation levels 
Di = 2Do that are assigned to the sites in a random 
fashion along the backbone so that the translational 
invariance of the lattice is broken by disorder. For the 
simulations presented by us in the manuscript, that 
were performed with the same realization of disorder 
for all considered temperatures, the breather moves 
towards the right end of the DNA chain. However we 
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Figure 5. The case T = 100 K: Spatio-temporal pattern of 
the inter-base distance dn(t). Parameters as in Fig. El For 
better illustration only a segment of the DNA lattice is 
shown. 
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Figure 6. Influence of temperature on the time-evolution of 
the breather center. Temperature as indicated on the curves. 



found that there exist also reahzations of the random 
Do — Di sequence for which the motion heads towards 
the left end. The fact that the breather movement is 
not diffusive, i.e., with random changes of directions, 
means that the moving breather is a stable new entity 
produced by the perturbation of the standing one. In 
mathematical terms, the thermal bath brings about 
the crossing of a separatrix between the basins of 
attraction of both breathers. Varying the temperature 
yields qualitatively equal results, i.e. the radial and 
torsional breather become eventually mobile traveling 
towards one end of the DNA lattice. Generally, it holds 
that the higher the temperature the more grow the 
amplitudes of the breather and the stronger the helix 
unwinds. For ambient temperature, i.e. T = 300K, the 
maximal amplitude of the radial breather is ten times 
larger than the one of the zero-temperature case being 
related to stretching of the associated hydrogen bond by 
0.6 A away from its equilibrium value. Moreover, at the 
sites (base pairs) where the helix is radially stretched 
to the most extend it is unwound by an angle of 12 °. 

Quantitatively, the results regarding the temperature 
dependence of the breather evolution are suitably sum- 
marized by the time-evolution of the first momentum of 
the energy distribution defined as 



N 



n{t) = ^ ^{uc-n) En^,{t) , 



(20) 
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and the energy En.i is defined in Eq. Q and ric is 
the site index corresponding to the center of the DNA 
lattice. 

This quantity describes the temporal behavior of the 
position of the center of a breather. Thus it represents 
a measure for the mobility of the breathers. Generally, 
the higher the temperature the larger the amplitude of 
the radial breather becomes and the faster the radial 
(and torsional) breather travels along the DNA lattice 
which is illustrated in Fig. As the degree of energy 
localization in the breathers is concerned, the energetic 
partition number serves as an appropriate mean. It is 
defined as 
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Figure 7. Temperature dependence of the temporal behavior 
of the partition number. Temperature as indicated in the 
graph. 



P quantifies how many sites are excited to contribute 
to the breather pattern. In the T = case the partition 
number changes only slightly and performs oscillation 
in an interval being rather close to its initial value for all 
the time. The relatively small growth of P(i) accounts 
for the phonons being emitted from the central region. 

For T > we find that Pit) rises with time. While up 
to times i < 80 the mean growth of P(t) proceeds lin- 
early, afterwards, when the DNA lattice system comes 
more and more closer to equilibrium, the participation 
number increases only in a logarithmic fashion. Even- 
tually, the value of P(t) reaches a plateau around the 
time value of 150. Compared to the initial width of the 
localized distortion pattern, i.e. P(0) = 29, the energy 
spread has risen by < 44%. This energy distribution 
over the lattice, however, is mainly due to the thermal- 
ization of the lattice by the thermal perturbation rather 
then by an actual spread of the radial breather. The 
latter maintains its degree of localization throughout 
the time. The time evolution of P is shown in Fig. 
Remarkably, in the range of 50 K ^ 300 K we observe 
almost equal temporal evolution of P regardless of the 
value of T. Thus the thermal perturbations do not 
strongly influence the width of the breather anymore 
giving evidence of the stability of the breather and the 
pronounced storing capability of the DNA lattice 
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4. Discussion 

In conclusion, we observed that out of an initial non- 
equilibrium situation, for which the hydrogen bonds of 
an under-twisted segment of the DNA lattice have been 
stretched, breathers develop in the radial and angu- 
lar displacement variables. For zero temperature the 
breathers remain standing in the initially excited region. 
Interestingly, for T > the breathers start to travel 
coherently towards an end of the DNA duplex. The 
observed breathers represent realistically the oscillating 
bubbles found prior to complete unzipping of DNA, i.e. 
they oscillate with periods in the range 0.3 — 0.8 ps, 
are on the average extended over 10 — 20 base pairs and 
possess maximal amplitudes of the order of < 0.6 A. Our 
results demonstrate that the action of some enzyme, 
mimicked by localized radial and torsional distortions 
of the DNA equilibrium configuration, initiates in fact 
the production of oscillating bubbles in DNA. Moreover, 
these oscillating bubbles sustain the impact of thermal 
perturbations. Whether the bubbles, when meeting each 
other, can merge such that a larger-amplitude radial 
breather results which resembles the replication bubble 
with broken H-bonds is still being investigated. 
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